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Abstract. - In this paper, we relate the coupling of Markov chains, at the basis of perfect sampling 
methods, with damage spreading, which captures the chaotic nature of stochastic dynamics. For 
two-dimensional spin glasses and hard spheres we point out that the obstacle to the application 
of perfect-sampling schemes is posed by damage spreading rather than by the survey problem of 
the entire configuration space. We find dynamical damage-spreading transitions deeply inside the 
paramagnetic and liquid phases, and we show that critical values of the transition temperatures 
and densities depend on the coupling scheme. We discuss our findings in the light of a classic 
proof that for arbitrary Monte Carlo algorithms damage spreading can be avoided through non- 
Markovian coupling schemes. 



Introduction. — Chaos manifests itself in Hamilto- 
nian dynamical systems when any two nearby initial con- 
figurations drift apart with time. Chaos can also be de- 
fined for cellular automata and for Markov chain algo- 
rithms. In these dynamical systems, following Kauffman 
PP, the drifting-apart of configurations is termed "dam- 
age spreading". In contrast, for "regular" dynamics, two 
nearby initial configurations become identical after a fi- 
nite time, and remain indistinguishable from then on. For 
Markov-chain Monte Carlo algorithms, the closely related 
case where the entire space of initial configurations be- 
comes identical is termed "coupling". Once it has cou- 
pled, the Markov chain has lost all correlations with the 
initial configuration. The coupling of Markov chains has 
risen to great prominence when Propp and Wilson used it 
for a perfect sampling method for Markov chains named 
"Coupling From The Past" (CFTP) [2]. When applicable, 
this method overcomes the problem of estimating the cor- 
relation time of a Monte Carlo calculation. In the present 
article, we shall discuss the fruitful connection between 
damage spreading and coupling [3]. 

In systems with N elements (spins, hard spheres, etc), 
the configuration space generally grows exponentially with 
N, and CFTP thus faces two distinct challenges. First, 
it must survey the entire configuration space in order to 
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prove coupling. Second, it must avoid damage spreading 
which would cause the coupling time to explode: it would 
become much larger than the correlation time as any two 
configurations have a very small probability for finding 
each other in a large space. 

The surveying problem is avoided in systems with a spe- 
cial property called "partial order" , as for example the fer- 
romagnetic Ising model under heat-bath dynamics JHH]. 
For more general systems (without partial order, but with 
local update algorithms), such as spin glasses and hard 
spheres, a recent "patch" algorithm inspired by numerical 
block scaling ideas OH] allows us to rigorously follow a 
superset of all initial conditions until it couples. This al- 
gorithm generates only modest overheads of memory and 
CPU time [3]. It was found that the coupling can be es- 
tablished after a time evolution very close to the coupling 
time. 

The second problem, the explosion of the coupling time 
related to damage spreading, poses the veritable obsta- 
cle to the application of CFTP ideas. Damage spread- 
ing has been studied in many physical systems, in par- 
ticular spin glasses j6]. In several spin glass models with 
heat-bath dynamics, it is now well established that a dy- 
namical damage-spreading transition occurs at a critical 
temperature, /?d s , located in the paramagnetic phase [7]: 
the dynamic is regular at temperatures higher than l//3ds 
and chaotic at lower temperatures. Even for the two- 
dimensional ±J Ising spin glass, which has a thermody- 
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namic phase transition at T = [SHTT]. the transition to 
chaos takes place at finite temperature [7J. We study in 
this paper the damage spreading of spin glasses and hard 
spheres, the divergence of the ratio of the coupling time 
and the correlation time, for different algorithms. 

Random walks in high dimensions. Before an- 
alyzing two-dimensional spin glasses and hard spheres, 
we illustrate coupling and damage spreading in a simple 
Markov chain algorithm that can be interpreted either as 
a random walk in an A-dimcnsional hypcrcubic lattice, 
as the dynamics of N distinguishable non-interacting par- 
ticles in a one-dimensional lattice of length L, or as N 
non-interacting Potts spins with L states. For the ran- 
dom walk (see Fig. [1} , each A-dimensional lattice site i = 
{io, . . . , iN-i} is described by integers ik € {0, . . . , L — 1} 
with periodic boundary conditions in the ik ■ The parti- 
cle can hop from site i to one of i's nearest neighbors in 
direction k, j = i ± 5k, with 5k = {0, ...,1,0,...} (pe- 
riodic boundary conditions are again understood). The 
probability for moving from i to j is 



p(i ->■ j) = 



i 

3N 
1 
3 



for j 
for j = i 
otherwise 



i ± 5 k 



(1) 



The simulation thus samples at each time step one dimen- 



time t, a randomly chosen particle k hops to the left or 
to the right, or it remains on the same site, each with 
probability 1/3, as above. 

A two-configuration coupling is a random process p(i — > 
—> j') for the joint evolution of two random walks 
such that integrating over one of them yields the original 
random walk of Eq. ([T]) for the other. After they meet, the 
two configurations evolve in the same way. The simplest 
choice for a coupling is the product ansatz, 



p(i -> j, i -> /) = < 



p(i 
p(i 




j)v{i' 
J) 



3 



if i ^ i' 

if i = i', j = j' 

otherwise 



(2) 

where the two random walks evolve independently from 
each other if they arc on different sites i and i', but stay 
together once they have met (J = j' if i = i'). To imple- 
ment this coupling for any number of configurations, one 
samples at each time step independent random moves at 
each site, so that particles on the same site experience the 
same randomness. In the above-mentioned representation 
of particles on the one-dimensional line, we consider the 
coupling of two ./V-particle systems, again described by 
Eq. ([2]), as the independent evolution, at time t, of the 
L N possible configurations of the system. Naturally, the 
coupling time scales as L N whereas the correlation time 
(in sweeps) behaves as L 2 . 



coup 




Fig. 1: Coupling of two random walks in a periodic N- 
dimensional hypercubic lattice of length L. After the time 
Tcoup, the two random walks evolve identically. The chaotic 
coupling of Eq. ([2]) is shown. For the regular coupling of 
Eq. (J2J), the displacement at time t is in the same direction 
k, and it is a function of ik only. 



sion, k, among the N available ones (it moves in "x" , or 
"y" or "z", etc). In dimension fc, it then hops with prob- 
abilities 1/3 each to the left or to the right, or remains on 
the same site. Equation ([T]) also describes N distinguish- 
able non-interacting particles on a one-dimensional lattice 
of length L, again with periodic boundary conditions: At 



100000 



10000 



1000 




10000 
. 1000 



• V i,:i: l')000 - 

1000 - 
100 
10 




1000 



t/N 



coup 



a logN+b 

chaotic ir- 
regular -* 



10 



100 



1000 



N 



Fig. 2: Chaotic and regular couplings for the random walk 
in a iV-dimensional hypercubic lattice of length L = 5 (see 
Eqs © and (J3J), respectively). The random process for a single 
random walk is defined by Eq. flj in both cases. Inset: Number 
of configurations vs. time as a function of the number iVinit of 
initial configurations, for N — 6, L — 5 for the chaotic coupling. 
The same realization of the coupling of Eq. Q is used for all 



An alternative coupling consists in sampling, at time 
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t, one dimension k common to all random walks. The 
two-configurations coupling scheme is then 



if i k + i' k 




tfik = i' k ,jk=j' k (3) 
otherwise 



so that two configurations i and j with i k = i' k will pre- 
serve this common coordinate (j k = j' k ). In the rep- 
resentation of N particles on a one-dimensional lattice, 
the same particle k is selected for each configuration, 
and for two different configurations, the particles labelled 
k stay together once they have met on the same site. 
The dynamics is then regular and the coupling time is 
T"coup/^V ~ a log N (sec Fig. [2]). The logarithmic behaviour 
is explained by the fact that particles move independently 
from each other, the coupling time for the entire system 
is thus the maximum of the N coupling times for each 
particle. 

In conclusion, we see that the same iV-dimensional ran- 
dom walk of Eq. ([1]), with a correlation time of order L 2 , 
allows two very different coupling, chaotic and regular. In 
spin glasses and hard spheres, these regimes are realized 
for the same coupling at different temperatures. 

Spin glass. — The random walk considered previously 
can also be considered as an L-statc Potts model at infinite 
temperature evolving under heat-bath dynamics. 

The product ansatz of Eq. © would correspond to the 
independent evolution of the spin configurations, and it is 
clearly chaotic. With the coupling of Eq. ((3j), in contrast, 
all spins evolve and couple independently at /3 = 0, and the 
global coupling time is again the maximum of the coupling 
times of the individual spins. The Monte Carlo dynamics 
is thus regular, and the diagram of Fig. [5] carries over to 
the general case with L > 2. At finite temperatures /?, the 
energy of a spin configuration is given by 



We first consider heat-bath dynamics, which consists in 
choosing one spin Sk and updating it with probabilities 



7T( Sk = ±1) 



1 



1 + cxp (T2ft fc /3) ' 



(4) 



where the field on site k is given by hk = X)z Jkisi- The 
coupling is defined by the use of the same random numbers 
for each configurations. 

For the two-dimensional ferromagnetic Ising model (all 
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1,L — 2), the dynamics remains regular at all 



temperatures. Below the Curie temperature, r c 
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very large, but so is the correlation time r CO rr, and the 
partial order implies that the complexity of t coup /t coii 
< 0(\ogN) [2]. The partial order is preserved in the 
disordered Ising model with ferromagnetic interactions 



Jij = Jji > 0, and in this model also, the theorem of 
Propp and Wilson guarantees that r coup is, up to loga- 
rithms, of the same order as r corr . 

Frustrated models, as for example spin glasses, do not 
exhibit partial order, and can thus undergo a damage- 
spreading transition. In the two-dimensional ±J Ising 
spin glass, the quenched random interactions satisfy Jij = 
Jji = ±1 with equal probability. Although this model 
is paramagnetic for all finite temperatures, Campbell and 
de Arcangelis [7] found a damage-spreading transition for 
the heat-bath algorithm at /3ds — 0.59. In previous work 
3,5, we succeeded in coupling large systems down to this 
temperature using the patch algorithm. We showed that 
the patch algorithm's upper bound on the coupling times 
agrees well with the lower bound obtained from a partial- 
coupling approach, where one checks coupling for a finite 
number iV; n it of random initial conditions rather than for 
the entire configuration space (iV; n ;t = 2^). As shown in 
the inset of Fig. [3j for one realization of the random pro- 
cess, the coupling time does not vary if N- ln i t > 10, and 
for iVi n it = 1000 it equals the coupling time for the entire 
configuration space. 
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Fig. 3: Disorder- averaged coupling time for the heat-bath algo- 
rithm of the two-dimensional ±J Ising spin glass. A dynamical 
phase transition is seen at the damage-spreading temperature 
Pds — 0.58. Inset: Saturation phenomenon for TV = 64 2 spins 
at j3 = 0.56. 

In the main graph of Fig. [3j we show the coupling 
time as a function of N at constant temperature. A dy- 
namical phase transition is seen at the damage-spreading 
temperature /?d s — 0.58. In the chaotic phase, T coup /N 
grows exponentially with N, but only logarithmically in 
the regular phase. The dynamical phase transition in this 
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model (without a spin glass phase at finite /3), although 
not mathematically proven, appears firmly established. It 
is without influence on single-particle properties. To il- 
lustrate this point, we verify that the correlation time 
Tcorr/N, computed with the autocorrelation function 

1 N 

^) = ]v"5>* (0)s * W) ' ( 5 ) 

i=0 

remains constant in the chaotic phase and only t coup /t coit 
diverges with TV — > oo. 

After the heat-bath algorithm, we now discuss the 
Metropolis algorithm, where individual spins are 
flipped with a probability depending on their local field. 
In the standard implementation, spin flips are accepted 
with a probability equal to 1 at infinite temperature. To 
allow coupling at any (3 we use 

2 

p(sk -> Sk) = - min(l,cxp(-2/3s fc /i)). (6) 

At each step the same spin k is updated for all copies of the 
system. For this dynamics, several coupling schemes can 
be set up. If the same random number 7 is used for each 
configuration, the coupling does not take place, as two op- 
posite configurations will always stay opposite. We adapt 
the regular coupling of Eq. ([3]) and use two independent 
random numbers, 71 for "up" spins and 72 for "down" 
spins. The coupling time has then the same qualitative 
behaviour as in Fig. [3j logarithmic at high temperatures 
and exponential at low temperatures, but with a critical 
temperature /?d s — 0.33. The Metropolis algorithm, with 
this coupling scheme, has thus a higher dynamical criti- 
cal temperature than the heat-bath algorithm, with which 
it shares all the qualitative features. This confirms that 
the dynamic damage-spreading transition is algorithm de- 
pendent. One may also choose the random numbers in 
the Metropolis algorithm using 7 for Sk = 1 and 1 — 7 for 
Sk = — 1 • This scheme correlates opposite spins better and 
the critical temperature is found to be /3ds — 0.52. This 
results shows that, like for the previous random walk, the 
same Markov-chain allows for qualitatively different cou- 
pling. 

Hard spheres. — After spin glasses, we now con- 
sider another key model in statistical physics, namely hard 
spheres. This model's Hamiltonian dynamics, realized in 
the event-driven molecular dynamics algorithm [Ud^], is 
chaotic for all densities and for all N [4][l3l[Il]. In this 
section, we will present several Monte Carlo algorithms 
for hard spheres, which allow for coupling of the entire 
configuration space. Two of the algorithms remain reg- 
ular below a finite critical packing fraction, r/^s, in the 
limit N — > 00. In the following discussion we are not con- 
cerned with algorithmic efficiency of the implementation, 
and only concentrate on the coupling properties. 

Birth- and- death algorithm. In the grand-canonical 
birth-and-death Monte Carlo algorithm, particles are 



placed inside the box at random positions x = (xk,yk) 
at rate A if no overlaps with previously placed disks are 
generated. The life time of each disk is sampled from an 
exponential distribution with rate 1. One realization of 
the algorithm is represented in the diagram of Fig. |4l The 
mean number (TV) of particles in this system is controlled 
by the activity A. 




Fig. 4: Grand-canonical birth-and-death algorithm for hard 
disks. Disk i appears at time ti, at position Xj = (xi,yi), and 
it disappears at time ti + n. In time, disk i describes a cylinder. 
Disks (cylinders) which are accepted, because they create no 
overlaps with earlier disks, are drawn in dark gray, while re- 
jected disks are drawn in light gray. The configuration space of 
this system is infinite, yet the possible configurations at time 
t are a subset of the finite set (of dark and light cylinders) 
produced from a horizontal cut in this diagram. 

This model's state space is infinite, but the survey of all 
possible initial conditions is nevertheless feasible [51IT5]. 
For any realization of the algorithm, the possible config- 
urations at time t are a subset of the finite set produced 
from a horizontal cut in the diagram of Fig. [4] The patch 
algorithm again yields sharp upper bounds for r coup [5]. 
Surprisingly, this algorithm for hard disks remains regular 
below a finite density 77ds in the limit N — > 00 [151116] , 

We again study damage spreading in this model by ap- 
plying the same Monte Carlo dynamics (same choice of 
Xj,tj,Ti) to iVi n it random hard-sphere initial conditions at 
time t — with life times sampled from an exponential 
distribution. The data shown in Fig. [5] again indicate a 
dynamical phase transition between the regular regime at 
packing fractions 77 < 77ds — 0.29 and the chaotic regime 
above rj^s- This density corresponds to the limiting den- 
sity found with the patch algorithm [5]. 
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Fig. 5: Coupling time of the birth-and-death algorithm of Fig. [4] 
for two-dimensional hard spheres (estimated with iVi n it = 100). 
The damage spreading transition occurs at a packing fraction 
Vds ^ 0.29. 



sity 77d s — 0.13 (see Fig. [5]), which is smaller than for the 
closely related birth-and-death algorithm. 

Spot algorithm. We finally study the coupling for a 
Markov-chain similar to the Metropolis algorithm: the 
spot algorithm. The Metropolis algorithm for TV hard 




/ — t — 1 Tcoup — 2 



Fig. 7: Spot algorithm for hard spheres: The randomly chosen 
spot position defines the attempted move of a disk inside the 
spot. The spot radius satisfies o" Bpo t < a, and at most one disk 
is moved at time t. An example with N = 1 and a spo t = <r is 
shown. 



''Labelled displacement" algorithm. A canonical ver- 
sion of the birth-and-death algorithm is the "labelled dis- 
placement" algorithm where, at times t = 0,1,2,..., a 
randomly chosen particle k is moved to a random posi- 
tion Xfc, if this move creates no overlaps. We see clear 
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Fig. 6: Coupling time for the labelled displacement algorithm. 
The dynamical transition to chaos occurs at a lower density 
(ryds — 0.13) than for the birth-and-death algorithm of Fig. [S] 

evidence of a dynamical phase transition at a critical den- 



spheres consists in moving, at time t, a particle k by a 
random vector 5 = (6 x ,5 y ). As the configuration space 
is continuous, the coupling probability is zero if one uses 
a naive coupling scheme. The following spot algorithm 
is more successful (although we will show its coupling is 
chaotic at all densities): at time t, it places a spot, a 
disk-shaped region with radius fJ spo t < er, at a random 
position x s . The spot contains at most one disk center, 
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Fig. 8: Coupling time r CO rr of the spot algorithm for TV hard 
disks (estimated with JVj n i t = 100) as a function of packing 
fraction r\. For all rj, the coupling time is exponential in N, 
and we conjecture the coupling to be chaotic. 

and the move consists in placing this disk at x s , if this 
creates no overlap with other particles (See Fig. [7]). The 
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spot algorithm satisfies detailed balance, and it generates 
the same moves as the Metropolis algorithm. Moreover, as 
illustrated in Fig. [Jj it succeeds in coupling. However, as 
shown in Fig. [5J the coupling time of the spot algorithm is 
an exponential function of N for all densities: the coupling 
is always chaotic. 

Conclusion. — In conclusion, we have in this paper 
studied the relationship between the coupling of Markov 
chains, which is of critical importance for the subject of 
perfect sampling, and damage spreading, which exposes 
the chaotic nature of the Monte Carlo dynamics. 

For the two-dimensional ± J Ising spin glass, which lacks 
an equilibrium phase transition at finite temperatures, we 
confirm the existence of a dynamical phase transition at 
/3ds — 0.58 [7J for the heat-bath algorithm. For lower tem- 
peratures the coupling time explodes. The Metropolis al- 
gorithm has the same damage spreading behaviour but 
with higher critical temperatures: /?ds — 0.33 or /?d s — 0.52 
for two simple coupling schemes. All damage-spreading 
transitions for this system are deeply inside the paramag- 
netic phase. 

For the two-dimensional hard-sphere system, we ana- 
lyzed three local Monte Carlo algorithms, the birth-and- 
death algorithm, inspired from Poisson point processes, its 
canonical version (the "labelled displacement" algorithm), 
and the spot algorithm, a straightforward adaptation of 
the Metropolis algorithm. The first algorithm shows a reg- 
ular regime only for packing densities below rj^s — 0.29, 
the coupling time was then of the same order of magnitude 
as the correlation time. The canonical version of the birth- 
and-death algorithm had a critical density of rjds — 0.13. 
These transition densities are again deeply in the liquid 
phase. 

Both for spin glasses and for hard spheres, the rigorous 
survey of the configuration space [3] remains feasible for 
all temperatures and densities. The application of per- 
fect sampling methods to these challenging problems is 
thus not so much limited by the surveying problem, as the 
patch algorithm allows to track the evolution of the entire 
configuration space, but more by damage spreading, the 
underlying chaotic nature of the Monte Carlo dynamics. 

In this context, it is of great interest that Griffeath [TTJ 
has constructed a coupling that always remains regular: 
It realizes the coupling at time t and at position Xt of two 
Markov chains that have started at time t = at configu- 
rations Xq and X' with the minimum of the probabilities 
to go from Xo or from X' t to X t . Griffeath's coupling is 
non-Markovian and very difficult to construct in practice, 
but it may point the way to couplings that remain regular 
at lower temperatures and higher densities than the naive 
couplings we discussed in this paper. 
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